Scattering optics of snow 
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Permanent snow and ice cover great portions of the Arctic and the Antarctic. It appears in winter 
months in northern parts of America, Asia, and Europe. Therefore snow is an important component of 
the hydrological cycle. Also, it is a main regulator of the seasonal variation of the planetary albedo. 
This seasonal change in albedo is determined largely by the snow cover. However, the presence of 
pollutants and the microstructure of snow (e.g., the size and shape of grains, which depend also on 
temperature and on the age of the snow) are also of importance in the variation of the snow’s spectral 
albedo. The snow’s spectral albedo and its bidirectional reflectance are studied theoretically. The 
albedo also determines the spectral absorptance of snow, which is of importance, e.g., in studies of the 
heating regime in snow. We investigate the influence of the nonspherical shape of grains and of 
close-packed effects on snow’s reflectance in the visible and the near-infrared regions of the electromag- 
netic spectrum. The rate of the spectral transition from highly reflective snow in the visible to almost 
totally absorbing black snow in the infrared is governed largely by the snow’s grain sizes and by the load 
of pollutants. Therefore both the characteristics of snow and its concentration of impurities can be 
monitored on a global scale by use of spectrometers and radiometers placed on orbiting satellites. 
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1. Introduction 


Snow covers a large portion of the terrestrial surface 
during the winter season. Therefore snow fields 
have potentially significant effects on the planetary 
albedo and climate. The main factors in the influ- 
ence of snow fields are snow cover and snow albedo r. 
Snow is highly reflective in the visible, where the 
solar spectrum’s maximum is located. Then albedo 
ris close to 1. However, the albedo decreases in the 
UV and IR spectral regions. The decrease in albedo 
is the most pronounced in the near-IR spectral re- 
gion, with values ofr close to zero at 1.5, 2.0, and 2.5 
um, where ice absorption maxima are located. 
Snow is essentially black then. 

The snow’s albedo controls the surface energy bud- 
get and the rate of melting. Therefore, knowledge of 
the albedo helps in understanding the seasonal 
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change in snow cover. A large amount of both ex- 
perimental and theoretical research is devoted to 
studies of snow’s optical properties. Recent compre- 
hensive reviews!? contain ~300 references to the 
subject. One of the most important results that was 
found is the fact that snow’s spectral albedo can be 
fitted by radiative transfer calculations under the 
assumption of spherical shape of snow grains3-5 (see, 
e.g., Fig. 10 of Ref. 4 and Fig. 9 of Ref. 5). This fact 
has allowed algorithms for retrieval of snow’s grain 
size from airborne and spaceborne optical instru- 
ments to be developed.16 

It is known, however, that snow in situ has an 
extremely complex structure and consists of irregu- 
larly shaped ice grains in contact with one another 
(see, e.g., Fig. 5 of Ref. 2). The optical properties of 
particles are controlled not only by the size but also 
by the shape of the particles.” Therefore it is clear 
that the model of spherical particles has a limited 
applicability to snow optics. An examination of this 
physical problem is the main task of this paper. 
Note that some important advances in this direction 
have already been reported.® In particular, it was 
shown? that the actual shape of grains has a profound 
effect on snow’s bidirectional reflectance in the visi- 
ble. Another issue that is addressed in this paper is 
the influence of close-packed effects!© on snow’s re- 
flective properties. 

According to our main ideas about snow optics that 


1 March 2004 / Vol. 43, No. 7 / APPLIED OPTICS 1589 


are developed in detail in this paper, the snow optics 
should be constructed with use of the following: 


e A model of snow as fractal close-packed ice 
grains instead of a noninteractive spheres as in the 
conventional spherical model, 

e Simple geometrical optics equations instead of 
Mie calculation for the snow’s local optical character- 
istics; and 

e Special exponential, analytical, asymptotic so- 
lutions of the radiative transfer theory for global 
snow’s optical properties instead of computations 
with radiative transfer codes. 


Implementation of these ideas provides simple an- 
alytical solutions and analytical comprehensive algo- 
rithms for retrieval of snow’s grain size and pollution 
concentration. Such an approach has been devel- 
oped and carefully tested but poorly publicized until 
now. The main ideas of this approach were dis- 
cussed in several workshops on the retrieval tech- 
niques for the Global Imager instrument developed 
by the National Space Development Agency of Japan 
in 1996, 1997, and 1998. 

In this paper we give a comprehensive pattern of 
this new approach to snow optics. The paper is 
organized as follows: In Section 2 we discuss snow’s 
global optical characteristics such as spherical and 
plane albedos. The approximate formula for snow’s 
bidirectional reflectance is also given there. Section 
3 is devoted to studies of snow’s local optical charac- 
teristics (e.g., transport and absorption path lengths). 
The problem of determining snow’s grain size from 
reflectance spectroscopy data is considered in Section 
4. 


2. Global Optical Characteristics of Snow 


A. General Equations 


Size d of snow grains is usually close to 1 mm.? This 
means that the methods of geometrical optics811-12 
can be applied to studies of snow’s local optical char- 
acteristics in the visible and the near IR. It is easy 
to show that the photon free-path length in snow, L, in 
the spectral range considered is close to the value of 
d: \¥<Il~d¥ 

The radiative transfer equation cannot be applied 
if à œ 1.8 Then the simple geometrical optics 
scheme that lies at the very heart of the radiative 
transfer theory would no longer be applicable. As 
we shall see, such is fortunately not the case for snow 
in the visible and the near IR. Therefore we assume 
that the radiative transfer equation can be used for 
calculating snow’s radiative characteristics. 

In particular, we use the following approximate 
analytical solution of this equation, which is valid in 
the limit of small absorption4-15; 


R(q) = Ro(q)exp| -af (q)]. (1) 


Here R is the bidirectional reflectance (or the reflec- 
tion function!) of a semi-infinite snow layer. 
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Vector—parameter q has coordinates ðo, ð, and ọ, 
which are the incident zenith angle, the observation 
zenith angle, and the relative azimuth, respectively. 

The function f (q) in Eq. (1) is given by the following 
ratio: 


K (Bo) Ko(9) 


a (2) 
Ro(Bo, ð, $) 


f(q) = 


where Ry equals the value of R at zero absorption 
and!2 


3 3 Qu a/2 
0 


0 


X Raða ©, e)cos? © sin Sd9. (3) 


The function Ko(ðo) is called the escape function in 
radiative transfer theory. It determines the angular 
distribution of light escaping from the semi-infinite 
nonabsorbing medium in the framework of the Milne 
problem (with sources located at infinity inside a me- 
dium!6), 

We also have for a in Eq. (1) (Ref. 12) 


a=4 (lea? Sl aa: (4) 


where l »s is the absorption path length and lą is the 
transport path length. Studies of the accuracy of 
Eq. (1) compared with the solution of the radiative 
transfer equation can be found in Refs. 17 and 18. 
Equation (1) can be applied with high accuracy for 
values of a less than unity. 

Note that the phase coherence of scattered light is 
lost over length scales that are longer than lą. The 
path of a scattered photon can be then described as a 
random walk with step size lą. On length scales 
that are long compared with /,,. it is useful to regard 
this random walk as photon diffusion with diffusion 
coefficient D = vi,,/3, where v is the speed of light in 
the medium. Therefore lą is an important parame- 
ter ofour theory. The value ofl,»s is the length scale 
at which absorption effects become important. This 
parameter is large in the visible range of the electro- 
magnetic spectrum. However, it rapidly decreases 
with wavelength owing to stronger light absorption 
by ice in the near-IR region. In fact, only these two 
parameters of a scattering medium enter into the 
light-diffusion equation.12 

We have from Eq. (1), as a — 0, 


R(q) = Ro(q) — aKo(8o) Kð), (5) 


which is a familiar result of the radiative transfer 
theory.1° Equation (5) allows us to establish the 
physical meaning of a. That is, by performing the 
integration with respect to So, 3, and » we obtain 
from Eq. (5) for the spherical albedo! 


r=l1-a, (6) 


where we have used Eq. (3) and the energy conser- 
vation law [see Eq. (A5) in Appendix A]. Recall that 
spherical albedo r is defined as!& 


p 2m ay T 
| ae | sae | ieee ao. 


0 0 0 


where we have introduced the cosines = cos ðo and 
m = |cos 9}. 

The value of r gives the fraction of light energy 
reflected from a semi-infinite medium under diffuse 
illumination. The value of a is, therefore, the frac- 
tion of the absorbed energy at the same conditions 
(and r close to unity). It follows [see Eq. (4)] that 
media with more absorbing particles (smaller laps) do 
not necessarily give the largest absorption by a me- 
dium as a whole. Instead, it is the ratio Z = lą/labs 
that is of importance. In particular, different semi- 
infinite weakly absorbing disperse media that have 
the same values of Z have values similar to reflec- 
tances R(q), provided that functions R,(q) are close 
as well. The physical reason for this fact is quite 
clear. Indeed, it is not only the absorption by a local 
volume that matters but also the type of scattering 
diagram. In particular, scattering diagrams that 
are highly extended in the forward direction with 
large values of lą lead also to a larger number of 
scattering and absorption events (hence, to a larger 
probability of absorption by a medium as a whole) 
before the escape of photons from the medium. 


B. Escape Function 


Let us consider escape function K,(8) [Eq. (8)] in 
more detail now. It follows from the general radia- 
tive transfer theory that the function K,(8) describes 
the angular distribution of light emerging from a 
semi-infinite nonabsorbing layer with light sources 
placed at an infinite depth inside the medium (the 
Milne problem).?°?!_ This means, in particular, that 
this function should not depend strongly on the type 
of scatterer. Such a property of the function K(9) is 
confirmed by numerical calculations.1921_ In partic- 
ular, the following approximation can be used for all 
types of scatterer??: 


K,(€) = (3/7) + 28) (8) 


if € = 0.2.1215.22, Kquation (8) is not influenced by a 
particular single-scattering law and can be used for 
media of various types, including snow. 

We emphasize that Eq. (8) also describes the an- 
gular distribution of light transmitted by thick layers 
of snow.!2 This suggests a simple way for experi- 
mental verification of Eq. (8) (e.g., performing angu- 
lar measurements of transmitted light in a cold 
chamber). Such data are not available to us. Sowe 
have chosen an indirect way to check Eq. (8) as ap- 
plied to snow optics, as follows. 

One can obtain from Eq. (1) for a plane albedo r,,(€) 
(Appendix A) 


r,(€) = exp[—oK,(€)], (9) 


where r,,(€) is defined as 


1 20 
rA(é) = | 


0 


1 
de | R(é, n, ¢)ndn. (10) 


0 


A similar equation can be derived for the spherical 
albedo given by Eq. (7) (Appendix A). It follows that 
r = exp(—a). 

The quantity r,(€) has been measured by many 


authors (see, e.g., Refs. 4, 23, and 24). It follows 
from Eq. (9) that 
1 1 
K((é) = — In —. (11) 
a ae lE) 


The precise value of a is usually unknown for a given 
experiment. However, a does not depend on é. 
Let us define the normalized escape function K,,(é) 


= K(é)/K,(é*). Then it follows from Eq. (8) that 
K pee (12) 
on(&) z 1 T 2¢* , 


where &* is a cosine of the selected incidence angle. 
We assume that é* = 0.7. Clearly, any other value 
of &* > 0.2 [see Eq. (8)] can be used for the normal- 


ization. It follows from Eq. (11) that 
In r,(&) 
K S Daag 1 


Therefore, measurements of r,,(€) can be used to es- 
tablish the relative angular dependence of the func- 
tion K,(é). The remaining proportionality constant 
can easily be found from the condition 


2 i K,(é)édé = 1, (14) 


0 


which follows from the energy conservation law.19 

A comparison of the function Ko,,(é) derived from 
Eq. (12) at &* = 0.7 and of measurements of r,,(€) [Eq. 
(13)] is given in Fig. 1. We can see that the experi- 
mental data closely follow the approximation in the 
rangeé=0.4. The difference between curves can be 
attributed in particular to the inaccuracy related to 
taking experimental data from plots r,,(€) (especially 
at € < 0.4, where r, is close to 1). We also show the 
results of the exact calculation?! of the function Ko,,(é) 
for media characterized by the Heney—Greenstein 
(HG) phase function p(0) = >7_, g’P,(cos 9) in Fig. 1. 
Here P,(cos 0) is a Legendre polynomial. It is as- 
sumed that g = 0.75. Note that p(0) gives the prob- 
ability of photon scattering in the direction specified 
by scattering angle 0 and that g = % fg p(6)sin 0 cos 
§ d9 is the asymmetry parameter. 

Aoki et al.5 found that the HG phase functions are 
more suitable than Mie phase functions for studies of 
snow’s optical properties. It should be stressed, 
however, that the function of K,(é) is highly robust to 
variations in phase function [at least for = 0.2 (Ref. 
19)]. The exact result and the approximation in Fig. 
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Fig. 1. Dependence of the normalized escape function on the co- 
sine of the incidence angle according to measurements in Refs. 23 
and 24, Eq. (12), and the exact calculations?! for the HG phase 
function with g = 0.75. 


1 are in close correspondence at é = 0.2 (the error is 
less than 2%). Therefore we conclude that Eq. (8) 
can be used in conjunction with Eqs. (1) and (2) for 
studying the influence of absorption on the bidirec- 
tional reflectance of snow. Note that K,(é) increases 
with €; this leads to a decrease ofr, with € [Kq. (9)], as 
has been confirmed by field measurements.*2*.24 

We can see, therefore, that the plane albedo of 
snow is larger for the Sun at the horizon (small values 
of £) than for the overhead Sun. The physics behind 
this effect is quite transparent. Indeed, photons in- 
jected into the medium along slant paths have more 
chances to survive than those at perpendicular inci- 
dence. 


C. Bidirectional Reflectance 


To model the bidirectional reflectance function of a 
plane-parallel semi-infinite snow layer R(q) we 
should also know function Ro(q) for nonabsorbing 
snow [Eqs. (1) and (2)]. That function depends on 
phase function p(@) only. Various assumptions 
about function p(0) (see, e.g., Refs. 25-31) lead to 
diverse results as far as function Ro(q) is concerned. 
This fact was clearly demonstrated in Ref. 9 (see Fig. 
5 of Ref. 9) for three shapes of particles, namely, 
spheres, hexagonal cylinders, and fractals (see also 
Fig. 2 ofthis paper). Therefore it is of importance to 
make the correct selection of the phase function in 
modeling the bidirectional reflectance of snow. 
Several theoretical phase function models are 
given in Fig. 3, together with experimental data?> for 
acrystalline disperse medium. The spherical model 
should be rejected. In particular, if this model were 
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Fig. 2. Reflection functions of a semi-infinite nonabsorbing me- 
dium with spherical and fractal grains at nadir illumination as a 
function of observation angle. Calculations were performed by 
use of the radiative transfer code described in Ref. 9 for a refractive 
index n equal to 1.31 and a wavelength equal to 0.55 um. The 
model of fractal particles that we used is described in Ref. 29. A 
size distribution of spherical particles is given by the cloud model 
C.1 with effective radius aef = 6 um (Ref. 12) and a coefficient of 
variance of the particle size distribution equal to 7/7. The in- 
crease in the size of the particles has an insignificant influence on 
the results.19 


valid, rainbows and glories should be observed (see 
the solid curve in Fig. 3 at scattering angles 6 close to 
145° and 180° and the solid curve in Fig. 2 for angles 
at 35° and 0°). Clearly, this model is not valid for 
snow. The hexagonal cylinder model studied in Ref. 
9 should also be rejected because such shapes do not 
appear in snow in large numbers.? 

The random particle model?®27 and the fractal- 
particle model?8:29 are good candidates for modeling 
bidirectional reflectance, however. The first model 
is based on the distortion of a spherical shape by the 
Gaussian distribution of its surface with various pa- 
rameters. The fractal particle is constructed in the 
following way: The initial shape is a regular tetra- 
hedron. A part of its surface is replaced by a re- 
duced version of the same tetrahedron. The 
repetition of this procedure at smaller triangles leads 
to a fractal surface with the dimension In 6/In 2 ~ 
2.58.29 Further details of the stochastic-particle 
models and the numerical procedures for finding p(0) 
can be found in Refs. 26-30. 

We stress that functions p(6) for fractals?*-29 and 
random particles?®?’ are similar (Fig. 3). Therefore 
functions R,(q) for fractals and for random particles 
do not differ much. This fact underlines the fact 
that what is of importance is not the particular sto- 
chastic particle model but rather the level of random- 
ness in a given model. 

We choose the fractal model here because it has no 


fractal model 

random particle model 
spherical particle model 
measurements 


phase function 


20 40 60 80 100 120 140 160 180 


scattering angle(deg) 


Fig. 3. Phase functions of ice particles for various theoretical 
models and measurements.”> Calculations were performed with 
the geometrical optical approximation at the wavelength 0.55 um, 
where absorption of light by ice is extremely weak. The fractal- 
particle model is described in Ref. 29. The random-particle model 
is described in Ref. 26. A standard deviation of o = 0.2 and 
correlation angle T = 5° were assumed for calculations according to 
the random-particle model.26 The sizes of the particles in fractal- 
and random-particle models do not influence the results given 
here. The phase function for spherical particles was calculated 
for the same spherical polydispersion as for Fig. 2 but for the 
effective radius a.f = 15 pm. 


free parameters for nonabsorbing ice particles and 
corresponds closely to the random-particle model?¢ at 
extreme values of its randomness parameters.2¢27 
The final decision, however, should be made when 
measurements of the snow’s phase function in situ 
are performed. 

It should be stressed that the poor correspondence 
of measurements and the theory for fractals near 90° 
and 170° (Fig. 3) could be due to the presence of small 
crystals or even of ice spheres in the experiment?® 
(see the photomicrograph in Ref. 25). The authors of 
Ref. 25 have been able to describe their experimental 
curve with high precision by using the unified theory 
of light scattering by ice crystals, assuming the fol- 
lowing shapes of crystals: irregular bullet rosettes 
and plates with rough surfaces. Clearly, the combi- 
nation of different shapes can be used to model ob- 
served phase functions of crystalline disperse media 
with a high accuracy.2!_ However, here the problem 
of the uniqueness arises. Therefore we prefer to use 
a model of a single randomly oriented fractal particle, 
which has no free parameters attached. 

A comparison of measurements*? of R(%) at o = 0° 
and calculations that use Eq. (1) with Eq. (8) taken 
into account and function Ry given in Fig. 2 for fractal 


~-©— Hapke model 


O approximation 
+ experiment 


reflection function 


0 20 40 60 80 
observation angle (deg) 


Fig. 4. Reflection function of snow obtained from measurements?32 
made according to our approximation [Eq. (1)] and Hapke’s model. 


particles is shown in Fig. 4. The value of a was not 
specified in the experiment.?2 So we used a as a 
fitting parameter at angle 30° (Fig. 4). 

We can see that our model is capable of describing 
accurately the experimental bidirectional reflectance 
of snow. It has no free parameters except the value 
of R at ð = 30°, in contrast with Hapke’s model*? as 
applied to snow. In particular, it is assumed in the 
framework of Hapke’s model33-34 that 


R= ern + s(g8)]o (0) + H(é)H(m) — 1}, 
(15) 


where wọ is the single-scattering albedo, p(8) is the 
phase function, and 
By 
1 + (1/h)tan(g/2) ” 
Here h and By, are free parameters of the theory and 


g is the asymmetry parameter. Function H(é) can 
be calculated approximately with? 


s(g) = (16) 


A= 
1+€y1— wo 


For Fig. 4 we used the following values of parameters 
in Eqs. (15) and (16): Bo = 1.0, wọ = h = 0.995, and 
the HG phase function p(§) with g = 0.449. 
Domingue et al.34 found that this set of parameters 
gives best fits of the measurements derived from 
Hapke’s model presented in Fig. 4. 

It follows that both our snow model and that which 
was derived according to Hapke’s theory explain the 


(17) 


1 March 2004 / Vol. 43, No. 7 / APPLIED OPTICS 1593 


almost constant value of the reflection function in a 
range of viewing angles smaller than 50° (near-nadir 
observations). This fact has also been confirmed by 
experiment.?2 However, our theory gives better 
agreement with measurements at the range of angles 
45°—80°. Note that the use of a more-complex model 
of the phase function in Eq. (15) fails to improve the 
accuracy of Hapke’s model to any great degree.*4 

Parameter a obtained from fitting measurement 
data®2 by Eq. (1) and that calculated from values of g 
and w, of Hapke’s model {a = 4[(1 —  )/3(1 — g)|'”” 
(Ref. 12)} are similar. In particular, it is equal to 0.2 
in the framework of our model and is only 10% larger 
(0.22) when it is derived from the fitting techniques 
based on Hapke’s equation. 


3. Local Optical Characteristics of Snow 


A. Transport Path Length 


Until now we have considered mostly the global op- 
tical characteristics of snow such as its bidirectional 
reflectance and albedo. Let us turn to the descrip- 
tion of snow’s local optical characteristics related to 
Eq. (1). In particular, we wish to establish the rela- 
tion of parameters l ps and J, to the microphysical 
characteristics of snow. This will also lead to the 
relationship between a [Eq. (4)] and grain size d. 

Let us start from lą. This parameter is defined by 
the following equation!2: 


= 1 
Sall = g) í 


where Sext is the extinction coefficient and g is the 
asymmetry parameter. For the large particles con- 
sidered here, o,,, and g are given as 


Oext = (Ca) + (Caz IN, (19) 


Cs BCS) 
(Ge lt Ga) 


where N is the number of grains ina unit volume. G 
and D signify the contributions" 12 of geometric op- 
tics (G) and diffraction (D) processes to extinction 
cross section C,,,, scattering cross section Cea, and 
asymmetry parameter g. The angle brackets mean 
averages with respect to the geometrical character- 
istics of snow grains (size, shape, and orientation). 

It follows from Eqs. (18)—(20) and the assumption 
that Cea ~ Cext which is valid for the weakly absorb- 
ing media considered here, that 


lir (18) 


(20) 


1 
~ NL = QIN C 


where Q = (Co.a79%)/ (Cca 4g"). Note that we have 
also used? g” ~ 1. 

Equation (21) shows that the transport path length 
in snow is defined by the geometrical part of the 
scattered field. This field is affected little by the 
close-packed effects (but is affected significantly by 
diffraction).1°.1233 Therefore we can neglect consid- 


li 


(21) 
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eration of the close-packed effects in the calculation of 
l,, because of their mutual cancellation. This is an 
important point. 

It is known! !2 that 


(Ceca) = (2)/4 (22) 


for nonspherical nonabsorbing convex particles in a 
random orientation, independently of refractive in- 
dex m of grains and their exact shape. Here (>) is 
the average surface area of particles. 

It follows from Eqs. (21) and (22) that 


4 
-NO - (g®))’ 


where we have assumed that Q = 1. Such is exactly 
the case for nonabsorbing large spherical particles. 
Then g@ does not depend on the size of particles, and 
(Coca? 8?) = (Coca Xe"). We assume that this equal- 
ity holds approximately for snow as a whole. 

Equation (23) can be written in a slightly different 
form, introducing the volumetric concentration of 
snow grains, C,,: 


Lin (23) 


C, = NV), (24) 
where (V) is the average volume of grains. Then we 
have 

2d 
= 2 
b= a- @ NG,’ om 
where 
(V) 
d=6-—. (26) 
(2) 


Equation (26) can be taken as a definition of effective 
grain size d for such a complex medium as snow. 
Note that for monodispersed spheres d is equal to 
their diameter. It follows for polydispersed media 
composed of spherical particles that d = 2a,,, where 
Aer is the effective radius, defined as the ratio of the 
third to the second moment of the particle size dis- 
tribution. 

The value of C, is close to 1/3 for snow.? This 
means that lẹ} ~ 2d/(1 — (g°)) >> \ because d > d 
and 0 < g? < 1 and also points to the possibility of the 
application of the radiative transfer equation to the 
problems of snow optics. 

It follows from Eq. (25) that the ratio M of transport 
lengths for media with spheres to that for media with 
nonspherical particles with the same values of C, and 
d is given by the following formula: 


=— eri, (27) 


where s and n here and below denote spherical and 
nonspherical scattering, respectively. Note that 
8„ ~ 1/2 for fractal particles considered here.35 
The value of the asymmetry parameter for large 
monodispersed spherical ice particles as a function of 


0 200 400 600 800 1000 


1200 1400 
size parameter 


Fig. 5. Dependence of asymmetry parameter g on size parameter 
x, obtained from Mie theory for monodispersed spherical particles 
that have refractive index 1.31 — 107i. 


size parameter x = md/) is presented in Fig. 5. We 
can see that g ~ 0.89 for large ice spheres. This 
means that!2 g,% = 2g — 1 ~ 0.78, which gives, when 
Eq. (27) is taken into account, M = 2.27. Therefore 
the transport path length in the medium composed of 
fractal particles is approximately two times smaller 
than in a medium with spherical ice particles that 
has the same values of C, and d. In particular, we 
have lą = 2.271,.", which has important conse- 
quences for optical particle sizing in snow. We con- 
sider this issue in more detail in Section 4 below. 


B. Absorption Path Length 


Let us study absorption path length labs now. It is 
defined as the inverse value of absorption coefficient 
Tabs: labs = lf Taps where Oabs — N(C avs) and (Cabs) 
is the average absorption cross section per particle. 
In principle, (C ps) can also depend on the number of 
particles N for close-packed media such as snow 
fields.36:37 However, in situ experiments?8 show 
that snow albedo [and therefore a, see Eq. (9)] does 
not depend on N. The value of a ~ Vlaps/Lir is in- 
dependent of N if (C,,,) is not a function of N. 
Therefore we neglect the effects of close-packed me- 
dia in the calculation of l ps and write 


abs 


= u (28) 
w OO 


where we have accounted for Eq. (24). 
Let us derive an approximate expression for the 
value of C,,,, for a particle of an arbitrary shape, as it 


is needed in snow optics. For this we use the defi- 


nition of Caps! 
k n" * 3 
Cabs = |E | E€ (r)E(r)E (r)d r, (29) 
u dy 


where k = 27/}, X is the wavelength, V is the volume 
of a particle, |E] is the absolute value of incident 
electric field vector Eo, £” = 2nx is the imaginary part 
of the relative dielectric permittivity of a particle, m 
= n — ix is the complex refractive index of a particle, 
and E(r) is the value of the electric field inside a 


particle at the point with radius vector r. We as- 
sume that 
E(r) = Ey > a,(r)explid,(x)], (30) 
s=0 


which corresponds to the representation of E(r) by a 
linear combination of simple waves with amplitudes 
a,E, and phases $,. It follows from Eq. (80) that 


|E|? = |E]? X, a? exp(—2¢,"), (31) 


s=0 
where ,” = Im(¢,) and we have neglected the phe- 
nomenon of interference, which is of no importance 
for snow grains with d > à. Then, using the fact 
that ,” is a small parameter for weakly absorbing 
grains, we obtain 
|E)? = BOr)| Eo)’, (32) 


where Bæ) = E£, a,7(r). 
It follows from Eqs. (29) and (32) that 


6. tee) | e")d’r, (33) 


vV 


where we have used the equality 


| y(r) F(r)dr = væ | F(r)dr. (34) 


v 


The value of rọ is generally not known. 
We neglect the internal inhomogeneity of grains. 
Then it follows from Eq. (33) that 


Cbs = nByV, (35) 


where y = 477x(A)/) is the absorption coefficient ofice. 
Also, we have 


(Cabs) = By\V), (36) 


where B = nf) and we have assumed that (BV) ~ 
(BV). One can obtain from Eqs. (28) and (36) 
1 
libs Sinayi 
ByC, 
Let us check the applicability of Eq. (35) for a spe- 


cific case of spherical particles. To do this we calcu- 
lated the ratio B = C,,,/yV with Mie theory for 


(37) 
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Fig. 6. Dependence of absorption enhancement parameter B on 
size parameter x obtained from Mie theory for monodispersed 
spherical particles that have refractive index 1.31 — 107i. 


monodispersed spheres with refractive index m = 
1.31 — 107i at à = 0.55 um as a function of size 
parameter x = 7d/d. The results are given in Fig. 6. 
We see that indeed B depends only weakly on x if we 
neglect the effects of interference, which are of no 
importance for snow optics because of the nonspheri- 
cal shape of grains and their polydispersity. The 
numerical data for the value of B scatter from 1.27 to 
1.31 for cases relevant to snow optics problems (x = 
500-1200). This gives a variability of less than 3%. 
We take the value B ~ 1.27 for estimations in this 
paper, which is similar to that used by Bohren®9 (B ~ 
1.26). This means [Eq. (32)] that B = B/n ~ 1.0 
(with errors smaller than 3%) for ice spheres in the 
visible, where n ~ 1.31. Also, we have on physical 
grounds B —> 1 as n — 1 [or Cy, — yV (Ref. 12)]. 
This means that B is highly robust against variations 
in the refractive index for the range of n relevant to 
snow optics problems. 

The values of B for spheroids and hexagonal cylin- 
ders with various aspect ratios were tabulated in 
Refs. 12 and 40. In particular, it was found that B = 
1.2-2.1 for spheroids and B = 2.2-2.5 for hexagonal 
cylinders if it was assumed that the ratio of their axes 
is in the range 0.5-2.0. The average of B for the 
interval 1.2—2.5 is 1.85, which is close to our estima- 
tion of B ~ 1.84 for fractal particles obtained from 
geometrical optics Monte Carlo calculations. The 
details of the code that we used to find B for fractals 
are given elsewhere.2® The close correspondence of 
B for fractals and a mixture of particles that have 
various shapes is not surprising when one takes into 
account the fact that chaotic scattering and absorp- 
tion by media with particles of diverse shapes should 
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approach those of a single, extremely irregular par- 
ticle such as a fractal. 

Therefore we propose to use the value B = 1.84 in 
snow optics instead of B = 1.27 as for spheres. This 
proposal, however, should be checked against care- 
fully planned experiments for natural snow. We 
note that the value of B may also depend on the age 
of the snow and on temperature. 

Our estimations and those given in Refs. 12 and 40 
allow us to conclude that large, weakly absorbing 
nonspherical ice particles generally have larger ab- 
sorption cross sections than spheres of the same vol- 
ume (B,, > B,). This means also that labs? > labs- 

In particular, we have for ice particles laps ~ 
1.45l aps” if the volume concentration of snow grains is 
fixed. Recall that /,.° = 2.27/,.”. This means that 
the ratio /,,/l,,,, which is of importance for reflec- 
tance problems [Eqs. (1) and (4)] is 1.57 times larger 
for spheres than for fractal particles with the same 
values of the product yd. 

As can be seen from Eqs. (1) and (4), this ratio can 
easily be determined from reflectance measurements. 
It is proportional to the size of the grains [Eqs. (25) 
and (37)|: 


lg 2Byd 
Labs 3(1 ~~ (g%)) l 


In particular, we have la/labs ~ 2.45yd for fractal 
particles and l,,/l,,, ~ 3.85 yd for spheres. It follows 
from Eq. (88) that d = $1,,./yla,<;, where $ = 3(1 — 
(g%))/2B. Let us suppose that the ratio 1,,/ la). is 
exactly known (the ratio /,,/l,,, is determined from 
the snow’s reflectance when the absorption coefficient 
of ice, y, is known with a high degree of accuracy). 
Then the error of the grain-size estimation that is due 
to the spherical shape assumption is given by 6 = 1 — 
b./a Taking into account that $, ~ 0.26 and $,„ ~ 
0.408, we obtain that 5 ~ 0.4, which allows us to 
conclude that the uncertainty in the determination of 
the grain size that is caused by shape effects could 
reach 40%. Therefore retrievals of snow grain size 
based on Mie theory may greatly underestimate the 
value of d. This finding is a major result of our 
study. Note that in situ measurements also confirm 
that an optically equivalent (in terms of spectral al- 
bedo fitting) diameter of the snow grain size is 
smaller than that measured in situ. This problem 
has been addressed, e.g., by Aoki et al.,> who sug- 
gested that the spherical approximation gives not the 
value of d defined in Eq. (26) but rather a dimension 
of a narrower portion of broken crystals or the branch 
width of dendrites. 

The value of a ~ Vi;,/l,,, is approximately 1.25 
times larger for spheres than for fractals. We con- 
clude, therefore, that although values of laps and ly 
differ considerably for spheres compared with frac- 
tals (see above), the difference is more moderate for 
the value of a, which determines the albedo of snow. 
Nevertheless, we stress that the use of the Mie model 
for the calculation of « with Eq. (4) leads to an over- 
estimation of a by 25% compared with the fractal 


(38) 


model. This leads to an overestimation of the snow’s 
heating (larger a) if the spherical-particle model is 
used and also gives the same error in the value of 1 — 
r as inr — 1 [Eq. (6)]. 


4. Determination of Snow Grain Size by Reflectance 
Spectroscopy 


Equations (2), (4), (8), (25), and (37) can be used for 
rapid yet accurate estimations of the spectral reflec- 
tance of snow for various observation geometries and 
snow grain sizes. These estimates can be used for 
several applications, including retrieval of the prop- 
erties of aerosols and clouds over snow fields. 

However, the most promising area for application 
of our semianalytical theory is in the determination 
of snow grain size by use of data from airborne and 
satellite-based radiometers and spectrometers. 

Snow grain size is an important parameter, which 
can be used, in particular, as an indicator of the age 
ofsnow. Changes in grain size help to identify snow 
dunes, melt areas, and regions of blue ice. They 
often indicate changes in the snowpack’s energy bal- 
ance. 

The rate of grain growth is exponentially propor- 
tional to temperature. This means that changes in 
grain size are indicators of thermodynamic processes 
in the snowpack. The estimation of grain size is of 
importance for calculating the absorption of radiation 
in snowpack. Note that increased absorption of so- 
lar light by snow leads to increased probabilities of 
melting snow and avalanches. 

To start with, let us find the analytical relationship 
between the snow’s bidirectional reflectance and the 
snow’s grain size. For this we substitute Eqs. (4), 
(25), and (37) into Eq. (1): 


R = R, exp(—bf yd), (39) 
where 
4 2B me 
p= = |] CE (40) 
3 la - | 


We conclude from Eq. (39) that the spectral depen- 
dence of In R can be presented as a linear function of 
X(A) = VyA). In particular, it follows that In Rà) = 
s + uX(\), where s = In Rọ and v = —bfVd. Our 
estimations show that the spectral variability of pa- 
rameters s and v can be neglected. 

Note that it follows [see Eqs. (1) and (39)] that 


a =byyd. 


This means that the spectral reflectance of snow is 
determined by the absorption of light on the length 
equal to grain diameter d (e.g., by the parameter c = 
yd). This is similar in some respects to the conven- 
tional transmittance spectroscopy. Then the value 
ofc is directly measured and yis obtained as the ratio 
c/d, where d is the length of the measured cell. 
Substituting the values B ~ 1.84 and (gf) ~ 0.5 
into Eq. (40), we obtain for media with fractal parti- 
cles b,, ~ 3.62. It follows that, for spheres with B ~ 


(41) 
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Fig. 7. Same as in Fig. 6, except for the value of b. 


1.27 and (gf) ~ 0.78, b, ~ 4.53. Figure 7 confirms 
this finding for spheres. We have, therefore, b,/b, ~ 
1.25, the result that has already been quoted above. 
Therefore we conclude that the value of b for spheres 
is approximately 25% too high than that of fractal 
particles. Therefore smaller grains are needed to 
produce the same value of R [Eq. (39)] for snow made 
from ice spheres as for fractal snow. 
It follows from Eq. (39) that 


1 R 
= —— In| — }. 
E yb’f oF | 


This physically based equation can be used as an 
alternative to techniques for retrieval of snow grain 
size based on fitting procedures (see, e.g., Ref. 6). 
Equation (42) can also be used for estimates of the 
influence of various factors on the retrieved value of 
d. In particular, we obtain from Eq. (42) that the 
uncertainty 5 in the value of y translates into the 
same value of uncertainty in the value of d. 

It follows from Eq. (42) under the assumption that 
Ro ~ 1.0 (Fig. 2) for both spheres and fractals that 
d,,/d, = (b,/b,)* ~ 1.6. Therefore we have d, ~ 
1.6d,, and currently used remote-sensing techniques 
may substantially underestimate the sizes of snow 
grains. Also, we have 6 = 1 — d,/d, ~ 0.4, as was 
outlined above. Note that the necessity to account 
for atmospheric scattering above snow (see, e.g., Ref. 
5) complicates the estimation of ò given above. 

It should be emphasized that all parameters in Eq. 
(42) (except Ro) are given by simple analytical formu- 
las specified above. This allows for the determina- 
tion of d from the measured value of R if Ro is 
obtained from precalculated lookup tables for fractal 
particles. Note that Ry is close to 1 for fractal par- 


(42) 
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ticles at nadir observation (Fig. 2) and nongrazing 
solar angles. Approximate equations for Rọ can be 
found elsewhere. !2:15.22 

The information on Rọ is not needed if the plane 
albedo of snow is measured. Then we have [Eq. (9)| 


C In%(r,) 
~ y(1 + 2 cos 9p)?’ 


where C = 49/96? ~ 0.42 and we have accounted for 
the fact that b ~ 3.6 for fractals. This allows for the 
immediate determination of d from r, if solar zenith 
angle ðo is known. It follows from Eq. (48) that 
uncertainty 6 in the value of y leads to the same value 
of uncertainty in the retrieved grain size d, as was 
already shown above. 

An even simpler equation can be derived for d if 
spherical albedo r of snow is measured. Indeed, it 
follows from Eq. (A8) of Appendix A in combination 
with Eq. (41) that 


r = exp(—b yd) (44) 
and, therefore, that 


(43) 


1 2 
d= yb? In’ r. (45) 
Interestingly, Eq. (44) is transformed into the well- 
known formula39.41 


r=1-byyd (46) 


for small values of light absorption in snow. Note 
that Bohren®® proposed the value b = 6.0, which is 
somewhat larger than the corresponding value for 
fractal particles (see above). It is also larger than 
values of b for asymptotically large ice spheres (Fig. 
7). Clearly, only experiments can determine the 
possible range of change of this parameter in natural 
snow. Note that we have from Eq. (46) that b = (1 — 
r)/Vyd for weakly absorbing media. This allows us 
to find the value of b from experimental measure- 
ments. 

The accuracy of Eq. (44) is studied in Fig. 8 by use 
of experimental data for snow reflectance obtained in 
Antarctica,*2 where snow pollution is of only minor 
importance. The data for refractive indices of ice 
were taken from Ref. 48. There is a problem with 
the choice of an appropriate value of d that is repre- 
sentative for measurements. Note that Grenfell et 
al.42 give vertical profiles of grain size measured vi- 
sually in the field at the South Pole (1985—1986) and 
the Vostok, Antarctica (1990—1991) stations. The 
average values of d were 0.15 mm for the South Pole 
station and 0.22 mm for the Vostok station (in the 
upper 10 cm of snow; see Table 3 of Ref. 42). They 
also determined average snow size profiles at the 
South Pole station (1986), obtained by the analysis of 
photographs (d = 0.28 mm in the upper 10 cm of 
snow), as presented in Table 4 of Ref. 42. The re- 
sultant average grain diameter was 0.22 mm, so we 
used this diameter in our calculations for Fig. 8. 
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Fig. 8. Dependence of spherical albedo on wavelength according 
to measurements?? and calculations made with Eq. (44) at b = 3.6 
and d = 0.22 mm. The data for refractive indices of ice are taken 
from Ref. 43. 


It follows from Fig. 8 that experimental data can be 
explained by the theory developed here for a wide 
spectral range starting from the UV until approxi- 
mately 1.4 um. For larger wavelengths the assump- 
tion of weakly absorbing media used extensively in 
this paper is not valid, so errors are increased. Then 
the exact radiative transfer equation with numerical 
calculations of the local optical characteristics of 
snow’s fractal grains should preferably be used. 
Such numerical results for a semi-infinite homoge- 
neous snow layer have already been reported.? The 
simplicity of the retrieval scheme as emphasized 
above is lost then, however. 

Clearly, such good fits as in Fig. 8 at) = 1.4 pm 
(and also for larger wavelengths) can be obtained 
with Mie theory as well. However, then the value of 
dis merely a fitting parameter. This parameter can 
differ substantially from the snow grain size defined 
by Eq. (26), as was emphasized above. This finding 
underlines the necessity of careful analysis of the 
spectral reflectance of snow relative to determination 
of the snow’s grain size. 

The analytical solution of the inverse problem pro- 
posed here can be applied to uniform snow only. 
Aged snow is characterized by the vertical inhomo- 
geneity of its physical properties, however. This 
means that the value of d derived from Eq. (42) will 
change with wavelength for vertically inhomoge- 
neous snow. This change can be advantageous, 
however. Indeed, in this case the reflection function 
is influenced by different parts of the layer, depend- 
ing on the wavelengths. Clearly, highly absorbing 
wavelengths will give information on the microstruc- 


ture of the top snow layer only. Deeper layers of 
snow can be sensed by weakly absorbing wave- 
lengths. This peculiarity was first studied by Li et 
al., who used data from the Airborne Visible/Infrared 
Imaging Spectrometer at wavelengths of 0.86, 1.05, 
1.24, and 1.73 um. Clearly, this issue deserves fur- 
ther exploration. For this exploration, however, our 
technique should be extended to vertically inhomoge- 
neous media as could be done, e.g., by use of Cham- 
berlain’s effective-reflective-layer technique.?!“4 
Another important issue that can be addressed in 
the framework of our approach is the identification of 
snow pollution (e.g., by soot). Then, assuming that 
the diameters of grains are obtained from IR chan- 
nels (see, e.g., Ref. 45), where y is determined almost 
exclusively by ice absorption, it follows from Eq. (42) 
that, in the visible, where ice is almost transparent, 


= 1 ml E 
YT Bfd aR) 


The value of y can deviate substantially from zero 
because of snow pollution. Clearly, the derived 
value of y can be used for estimation of the concen- 
tration of pollutants.45 Then the information on the 
type of pollutant (e.g., snow, dust, or both) should be 
known a priori or assessed from spectral curves R(\) 
[and related spectra y(A)]. 

Equation (47) offers a simple way to estimate the 
relative absorption coefficient of ice in snow y,(A) = 
y(\)/y(Ao), where ño is the selected reference wave- 
length. That is, it follows from Eq. (47) that 


(47) 


_ n{ERQ)] 


= DIER)’ (48) 


yA) 


where we have accounted for the fact that Z = 
1/Ro(ào) does not vary considerably with wavelength. 

An even simpler result follows from Eqs. (43) and 
(45): 


In? rA) Ièr) 
In? r Ao) In? rO) 


yA) (49) 


This allows for the direct determination of the rela- 
tive ice absorption coefficient from measurements of 
spectral albedo. 


5. Conclusions 


The science of snow optics is still in its infancy. Sim- 
plified models based on sparsely distributed ideal 
spherical particles are used studying the optical prop- 
erties of snow and also for deriving characteristics of 
snow by use of data from airborne and spaceborne 
spectrometers and radiometers.46-48 Such models 
are remote from reality. Snow is an extremely com- 
plex medium composed of nonspherical, irregularly 
shaped grains in contact with one another. This in- 
fluences the radiative and polarization characteris- 
tics of snow considerably. Nevertheless, it was 
found in previous studies*> that the spherical- 
particle model fits the measured spectral albedo re- 


flectance spectra of snow quite well. In this paper 
we have clarified the physics behind this result and 
shown that it cannot be considered a proof of the 
spherical model without independent control of the 
snow grain’s size. Also, the spherical approximation 
fails to explain the behavior of the snow’s bidirec- 
tional reflectance.5:32:34 We draw the reader's atten- 
tion to the fact that, unlike most conventional 
approaches, which compute the reflectance of snow 
with radiative transfer codes, this paper offers the 
use of an exponential asymptotic analytical equation 
for the reflectance of snow. This simple solution is 
shown to be highly accurate in the case studied. 

The main objective of this paper is to introduce and 
discuss a new approach to snow optics with a more- 
realistic model of snow as a medium with nonspheri- 
cal and close-packed snow grains. However, the 
results obtained can also be used far beyond the 
rather narrow topic considered here. In particular, 
they can be applied to other weakly absorbing light- 
scattering media with a complex structures, includ- 
ing whitecaps and bright soils. The last issue of 
particular importance for planetary spectroscopy.®? 

We propose to use a fractal-particle model instead 
of the spherical-particle model usually applied in 
snow optics. Clearly the fractal-particle model of 
snow is closer to reality than the spherical-grain 
assumption. We emphasize that the use of this 
more-sophisticated model does not lead to any com- 
plications (at least, in the UV, visible, and a consid- 
erable portion of the near-IR spectral regions, where 
absorption of light by snow is weak). The difference 
between these two models lies in the final equation 
for snow’s reflectance as different values of a form 
factor [b in Eq. (42)], which account for the shape 
effects. The value of this form factor could serve as 
a criterion with which to estimate and choose the 
appropriate snow model; it could and should be con- 
trolled in special experiments with natural snow. 

The use of the fractal model instead of the spherical 
model leads to ~40% larger values of the retrieved 
snow grain size compared with output from the 
spherical model. This means that retrieval tech- 
niques based on Mie theory will bias grain size [de- 
fined by Eq. (26)] considerably. This fact was 
underlined long ago by Grenfell et al.,49 who found 
disagreement (by 50%, which is close to our estima- 
tion) between measured and retrieved (in the frame- 
work of the spherical approximation) values of d. 
Similar results were reported in Refs. 5, 50, and 51. 
We believe that we have managed to establish the 
physics behind these findings. 

Note that the comparison of remotely sensed snow 
grain sizes with those determined in situ made in 
Refs. 46—48 also showed that remote-sensing tech- 
niques for measuring snow based on the spherical 
approximation do underestimate snow grain charac- 
teristics. It was argued that this underestimation 
could be due to snow pollution or to subpixel scene 
heterogeneity*® (e.g., unaccounted for contributions 
from thin snow layers, vegetation, soil, and rock). 
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However, the nonspherical shapes of snow crystals 
can also contribute to the discrepancy found. 

To avoid the influence of pollution (e.g., by soot) 
and atmospheric molecular absorption, one should 
measure snow reflectance in the near-IR spectral 
range outside gaseous absorption bands. Measur- 
ing in this range also minimizes the contribution of 
atmospheric scattering. Then one can use the ob- 
tained value of the snow grain size to find the level of 
light absorption in visible. The difference between 
this value and that of the absorption of pure snow 
gives the strength of snow pollution.4® Similar tech- 
niques can be used to identify various life forms in 
snow (see, e.g., Ref. 52). This could be of importance 
in the analysis of data from satellites orbiting other 
planets including Mars. Note that the unconfirmed 
possibility of life forms outside our terrestrial envi- 
ronment is a subject of intensive research.®?-54 

The analytical equations that we obtained can be 
applied for the solution of the inverse problem of the 
snow grain determination by space-based, airborne, 
and ground measurements. They can be also used 
to improve the snow hydrology for global circulation 
models.®® 

For the sake of simplicity, only horizontally homo- 
geneous snow layers of an infinite depth have been 
considered in this paper. The extension of this ap- 
proach to finite and to vertically inhomogeneous snow 
layers can easily be made by use of an approach 
similar to that described above (see, e.g., Refs. 12, 15, 
and 21). Remote sensing of heterogeneous snow 
fields (e.g., with soil and vegetation included in the 
scene studied) and issues related to the determina- 
tion of snow cover can be explored by use of a combi- 
nation of our approach with the linear or nonlinear 
multiple-end-member spectral mixtures analysis pro- 
posed in Ref. 48. 


Appendix A. Plane Albedo 
Here we derive a simple analytical equation for the 


plane albedo of a semi-infinite snow layer. The 
value of r, is defined by the following integral: 
1 2T 1 
lp = =| ae | ndyR(é, n, p), (A1) 
0 0 


where R(é, n, ¢) is the bidirectional reflectance, £ and 
y are cosines of the incident and viewing zenith an- 
gles, respectively, and ọ is the azimuthal angle. Let 
us substitute Eq. (1) into Eq. (A1). Then it follows 
that 


1 Qa 1 
fp > de ndnko(é, N, ẹ) 
T 


0 0 


x exp] = 


The value ofr, can be found from Eq. (A2) only by use 
of numerical integration. 


Pao Ka (A2) 


i RoE, N, ẹ) 
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Let us obtain an approximate analytical result for 
the plane albedo, using the expansion 


= | i ee _ yy ODK E)K) (A3) 
j RoE, 1, p) so s!R lé, n, ẹ) 
Then it follows from Eqs. (A2) and (A3) that 
r, = 1- aKo(¢), (A4) 


where we have neglected nonlinear terms in Eq. (A3) 
and accounted for the identities162° 


1 2a 1 
a dọ | ndynRo(é, N, ọ) a 1, 


0 0 


1 
2 | Ko(n)yndy = 1. (A5) 


0 


Unfortunately, our derivation cannot be extended to 
nonlinear terms in Eq. (A3). Therefore we assume 
that r, is governed by a similar law, as specified by 
Eq. (1): 


r, = exp[—aK,(é]. (A6) 


This formula is transformed into Eq. (A4) for small 
values of a. 

Equation (A6) can be used to find the spherical 
albedo, defined as [Eq. (7)] 


r=2 i r(é)édé. (A7) 


0 
Using Eq. (A6) and the same chain of arguments as 
above, we obtain from Eq. (A7) 
(A8) 


r = exp(—a). 


The accuracy of the derived equations was studied in 
Refs. 12, 15, 17, and 18. 
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